function [vox_coords vox_coords_times fig_handle fig_name fig_title] = oat_plot_roi_stats(S)

% [vox_coords vox_coords_times fig_handle fig_name] = oat_plot_glm_stats()
%
% Plot the copes and t-stats at a specified vox coordinate
%
% S.stats; % oat.first_level or oat.group_level results struct or fname
% (for which you will also need to specify S.oat)
% S.vox_coord; % vox index into oat stats OR string specifying operation on
% vox dimension, e.g. 'max_tstat'
% S.first_level_cons_todo; % list of first level cons to plot
% S.group_level_cons_todo; % list of group level cons to plot (if group
% level stats)
%
% MWW 2013

try, vox_coord=S.vox_coord; catch error('S.vox_coord not specified'); end;
try, stats=S.stats; catch error('S.stats not specified'); end;

if isstr(stats)
    stats=osl_load_oat_results(S.oat,stats);
else
    stats=stats;
end;

try, first_level_cons_todo=S.first_level_cons_todo; catch first_level_cons_todo=1:size(stats.cope,3); end;
try, group_level_cons_todo=S.group_level_cons_todo; catch group_level_cons_todo=1:size(stats.cope,5); end;

if ~isstr(vox_coord),
    if(isscalar(vox_coord))
        vox_coords=vox_coord*ones(length(first_level_cons_todo),length(group_level_cons_todo));
    else
        vox_coords=vox_coord    
    end;
else
    
    vox_coords=zeros(length(group_level_cons_todo),1);
    
end;

cols={'r','g','b','y','m','k'};
fig_handle=[];
fig_name={};
fig_title={};

for gconi=1:length(group_level_cons_todo),
    gcon=group_level_cons_todo(gconi);

    fig_handle(gconi)=sfigure;
    set(fig_handle(gconi),'Position',[1 1 1300 450]);
                    
    if(stats.level==2),
        fig_name{gconi}=['stats_tc_gcon' num2str(gcon) '_' stats.group_level_contrast_name{gcon}];
        fig_title{gconi}=['Time course of stats for gcon=' num2str(gcon) ', ' stats.group_level_contrast_name{gcon}];
        
    else
        fig_name{gconi}='stats_tc';
        fig_title{gconi}=['Time course of stats'];
    end;
    
    leg={};

    max_cons=6;
    if(length(first_level_cons_todo)>max_cons)
        disp(['Only plotting the first ' num2str(max_cons) ' contrasts']);
        first_level_cons_todo=first_level_cons_todo(1:max_cons);
    end;
    
    if isstr(vox_coord),
        % use first contrast to choose the voxel and tpts of interest
        con=first_level_cons_todo(1);

        cope=permute(stats.cope(:,:,con,:,gcon),[1 2 4 3 5]);
        stdcope=permute(stats.stdcope(:,:,con,:,gcon),[1 2 4 3 5]);
        tstats=cope./stdcope; %nvox x ntime x nfreq
        tstats=reshape(tstats,[size(tstats,1),prod(size(tstats))/size(tstats,1)]); %nvox x (ntime x nfreq)

        switch vox_coord                
            case 'max_abs_tstat',
                tstats=abs(tstats);
                [vox_coords(gconi) vox_coords_times(gconi)]=max2d(tstats);
                fig_title{gconi}=[fig_title{gconi} ' at voxel with max(abs(t-stat)) for ' stats.first_level_contrast_name{con}];
            case 'max_tstat',
                [vox_coords(gconi) vox_coords_times(gconi)]=max2d(tstats);             

                fig_title{gconi}=[fig_title{gconi} ' at voxel with max(t-stat) for ' stats.first_level_contrast_name{con}];
            otherwise
                error('Unknown vox_coord method');
        end;

    end;

    for coni=1:length(first_level_cons_todo),
        con=first_level_cons_todo(coni);
    
        cope=permute(stats.cope(vox_coords(gconi),:,con,:,gcon),[2 1 3 4 5]);
        stdcope=permute(stats.stdcope(vox_coords(gconi),:,con,:,gcon),[2 1 3 4 5]);
        
        % plot copes
        subplot(1,2,1);
        hold on;
        errorbar(stats.times,cope,stdcope,cols{coni});

        % plot 1-tailed tstats
        subplot(1,2,2);
        hold on;
        plot(stats.times, squeeze((stats.cope(vox_coords(gconi),:,con,:,gcon))./stats.stdcope(vox_coords(gconi),:,con,:,gcon)),cols{coni},'LineWidth',2); hold on;   

        leg{length(leg)+1}=stats.first_level_contrast_name{con};
    end;
    
    subplot(1,2,1);
    a=axis; axis([a(1) a(2) a(3) a(4)]);    
    plot4paper('time (secs)','cope'); 
    legend(leg,'Location','NorthWest'); 

    subplot(1,2,2);
    a=axis; axis([a(1) a(2) a(3) a(4)]);    
    plot4paper('time (secs)','t-stat'); 
    legend(leg,'Location','NorthWest'); 

    end;
end